#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy,math,scipy.stats
import psyco
from . import PCA,LSI
psyco.full()

def EuclidDistance(vecs):
    '''
    ユークリッド距離
    vecsはリスト（orタプル）のリスト。
    [[1,2,3],[2,3,4],[3,4,5],[4,5,6]]
    '''
    ret=dict()
    for i in xrange(len(vecs)):
        ret[i,i]=0
    
    for i1 in xrange(len(vecs)-1):
        v1=vecs[i1]
        for i2 in xrange(i1+1,len(vecs)):
            #すでに計算済みなら何もしない
            if(ret.has_key((i1,i2))):
                continue
            
            d=math.sqrt(sum( ((vecs[i2][i]-v1[i])**2 for i in xrange(len(v1))) ))
            
            ret[i1,i2]=d
            ret[i2,i1]=d
    return ret
    
    
def MahalanobisDistance(vecs):
    '''
    マハラノビス距離
    vecsはリスト（orタプル）のリスト。
    [[1,2,3],[2,3,4],[3,4,5],[4,5,6]]
    '''
    covinv=numpy.linalg.pinv(numpy.cov(numpy.matrix(vecs).T))
    matVecs=[numpy.matrix(vec) for vec in vecs]
    
    ret=dict()
    for i in xrange(len(vecs)):
        ret[i,i]=0
    
    for i1 in xrange(len(vecs)-1):
        v1=matVecs[i1]
        for i2 in xrange(i1+1,len(vecs)):
            #すでに計算済みなら何もしない
            if(ret.has_key((i1,i2))):
                continue
            dist=matVecs[i2]-v1
            #try:
            #    dist=matVecs[i2]-v1
            #except:
            #    print matVecs[i2]
            #    print v1
            #    import sys
            #    sys.exit(-1)
            d=math.sqrt((dist*covinv*dist.T)[0,0])
            
            ret[i1,i2]=d
            ret[i2,i1]=d
            
            del dist
    return ret

def PageRank(matSim,randomWalk=0.0,th=0.00001,x0=None,maxRepeat=1000):
    '''
    PageRank
    cite: S.Brin and L.Page: “The anatomy of large-scale hypertextual Web search engine”, Proc. of the 7th International World Wide Web Conf. 1998, pp.107-117 ( 1998)

    example:

    from yi01lib.Analyzer.LinkAnalyzer import PageRank
    print PageRank([[0,1,1,1,1,0,1]
                    ,[1,0,0,0,0,0,0]
                    ,[1,1,0,0,0,0,0]
                    ,[0,1,1,0,1,0,0]
                    ,[1,0,1,1,0,1,0]
                    ,[1,0,0,0,1,0,0]
                    ,[0,0,0,0,1,0,0]],0.1)
    #=>
    #[[ 0.28801181]
     #[ 0.16104133]
     #[ 0.13941992]
     #[ 0.1072456 ]
     #[ 0.18275003]
     #[ 0.05540376]
     #[ 0.06612755]]

    '''
    def normalize(lst,n=2):
        if isinstance(lst,list):
            l=math.pow(sum([v**n for v in lst]),1.0/n)
            return [v/l for v in lst]
        else:
            l=numpy.linalg.norm(lst,n)
            return lst/l
    floatall=numpy.vectorize(float)
    
    n=len(matSim)
    x0=[1.0/n]*n if x0 is None else normalize(x0,1)#初期値
    randvec=(numpy.matrix(x0)*randomWalk).T
    x=numpy.matrix(x0).T
            
    #A_T=numpy.matrix([normalize(lnk,1) for lnk in matSim]).T
    A_T=numpy.matrix(matSim).T
    A_T/=floatall(sum(A_T))
    
    for i in xrange(maxRepeat):
        x2=A_T*x
        if randomWalk>0:
            x2*=(1-randomWalk)
            x2+=randvec
        x2=normalize(x2,1)

        if  numpy.linalg.norm(x2-x)<th:
            break
        
        x=x2
    else: raise Exception("%d回の計算では無理" % maxRepeat)
    
    if isinstance(matSim,list):
        x=x.T.tolist()[0]
    return x

def HITS(matSim,th=0.00001,maxRepeat=1000):
    '''
    HITS
    cite: J. Kleinberg. Authoritative sources in a hyperlinked environment. Proc. 9th ACM-SIAM Symposium on Discrete Algorithms, 1998. Extended version in Journal of the ACM 46(1999). Also appears as IBM Research Report RJ 10076, May 1997
    example
    
    from yi01lib.Analyzer.LinkAnalyzer import HITS
    print HITS([[0,1,1,1,1,0,1]
                ,[1,0,0,0,0,0,0]
                ,[1,1,0,0,0,0,0]
                ,[0,1,1,0,1,0,0]
                ,[1,0,1,1,0,1,0]
                ,[1,0,0,0,1,0,0]
                ,[0,0,0,0,1,0,0]])
    #hub=[
    #[ 0.64642469],
    #[ 0.11208832],
    #[ 0.25505553],
    #[ 0.46620772],
    #[ 0.43118474],
    #[ 0.27395032],
    #[ 0.16186199]]

    #authority=[
    #[ 0.34668525],
    #[ 0.44219257],
    #[ 0.4991381 ],
    #[ 0.34840691],
    #[ 0.50063347],
    #[ 0.13940904],
    #[ 0.20899787]]
    '''
    n=len(matSim)
    x0=[1.0/n]*n#初期値    
    x=numpy.matrix(x0).T
    tA=numpy.matrix(matSim)
    A=tA.T
    def normalize(lst,n=2):
        if isinstance(lst,list):
            l=math.pow(sum([v**n for v in lst]),1.0/n)
            return [v/l for v in lst]
        else:
            l=numpy.linalg.norm(lst,n)
            return lst/l

    for i in xrange(maxRepeat):
        x2=normalize(A*tA*x)

        if  numpy.linalg.norm(x2-x)<th:
            break
        x=x2
    else: raise Exception("%d回の計算では無理" % maxRepeat)

    a=x2
    h=normalize(tA*x2)
    if isinstance(matSim,list):
        h=h.T.tolist()[0]
        a=a.T.tolist()[0]

    return h,a

def HITS4Community(matSim,th=0.00001,maxRepeat=1000,forceReturn=False):
    '''
    HITS
    cite: J. Kleinberg. Authoritative sources in a hyperlinked environment. Proc. 9th ACM-SIAM Symposium on Discrete Algorithms, 1998. Extended version in Journal of the ACM 46(1999). Also appears as IBM Research Report RJ 10076, May 1997
    example
    
    from yi01lib.Analyzer.LinkAnalyzer import HITS
    #4people, 7objects
    print HITS4Community([[0,1,2,1,1,0,0]
                ,[1,0,0,0,3,0,0]
                ,[1,1,0,0,0,0,0]
                ,[0,0,0,0,1,0,0]])
    '''
    def _normalizeRows(a):
        #縦に正規化
        return a/(1.0*sum(a))
    def _normalizeCols(a):
        #横に正規化
        return a/(1.0*a.sum(1))
 
    n=len(matSim)
    x0=[1.0/n]*n#初期値    
    x=numpy.matrix(x0).T
    A=numpy.matrix(matSim)

    step=_normalizeRows(A)*_normalizeCols(A).T

    for i in xrange(maxRepeat):
        x2=step*x

        if  numpy.linalg.norm(x2-x)<th:
            break
        x=x2
    else:
        if not forceReturn:
            raise Exception("%d回の計算では無理" % maxRepeat)

    h=x2
    a=_normalizeRows(A).T*x2
    return h,a


def SupportRank(A):
    '''
    A:support matrix
    =>どの人が(タテ)何をどれだけ(ヨコ)Supportしているか
    
    A=[[1,0,1],
       [1,2,1],
       [0,0,1],
       [0,1,1]]
    rankPerson,rankCluster=SupportRank(A)
    
    '''    
    zero1=numpy.vectorize(lambda o:0.000000000000000000000001 if o==0 else o)
    floatAll=numpy.vectorize(lambda x: float(x))    
    def kai2(O,E):
        '''カイ二乗検定量'''
        return numpy.multiply(O-E,O-E)/numpy.abs(zero1(E))
    _A=LSI(numpy.matrix(A),0.95) #categoryを軸に緩和
    _A=LSI(_A.T,0.95).T #Personを軸に緩和
    E=numpy.multiply(_A.sum(0),_A.sum(1))/zero1(_A.sum())
    lenPerson,lenPic=numpy.shape(_A)
    
    #0次picRank:単純に支持層の広さ
    
    #category
    #print scipy.stats.chisqprob(sum(kai2(_A,E)),lenPerson-1)
    goodPic=scipy.stats.chi2.ppf(0.95,lenPerson-1)/zero1(sum(kai2(_A,E)))
    badPic=1/zero1(goodPic)
    
    ##person
    ##print scipy.stats.chisqprob(kai2(_A,E).sum(1),lenPic-1)
    #print scipy.stats.chi2.ppf(0.95,lenPic-1)/zero1(kai2(_A,E).sum(1))
    O=_A/_A.sum(1)
    E1=goodPic/goodPic.sum(1)
    E2=badPic/badPic.sum(1)
    
    goodPerson=kai2(O,E2).sum(1)/zero1(kai2(O,E1).sum(1))
    badPerson=1/zero1(goodPerson)
    
    O=_A/sum(_A)

    E1=goodPerson/sum(goodPerson)
    E2=badPerson/sum(badPerson)

    goodPic=sum(kai2(O,E2))/zero1(sum(kai2(O,E1)))
    
    pictureRank=scipy.stats.f.cdf(numpy.sqrt(goodPic),lenPerson-1,lenPerson-1)
    personRank=scipy.stats.f.cdf(numpy.sqrt(goodPerson),lenPic-1,lenPic-1)
    personRank=personRank.T
    if isinstance(A,list):
        pictureRank=pictureRank.tolist()[0]
        personRank=personRank.tolist()[0]
    return personRank,pictureRank

def bimyoSupportRank(A):
    '''
    A:support matrix
    =>どの人が(タテ)何をどれだけ(ヨコ)Supportしているか
    
    A=[[1,0,1],
       [1,2,1],
       [0,0,1],
       [0,1,1]]
    rankPerson,rankCluster=SupportRank(A)
    
    '''    
    
    zero1=numpy.vectorize(lambda o:0.000000000000000000000001 if o==0 else o)
    floatAll=numpy.vectorize(lambda x: float(x))    
    def kai2(O,E):
        '''カイ二乗検定量'''
        return numpy.multiply(O-E,O-E)/numpy.abs(zero1(E))

    def rank(A):
        _A=PCA(A/floatAll(A.sum(1)),0.999)[1]
        E=numpy.multiply(_A.sum(0),_A.sum(1))/zero1(_A.sum())
        return scipy.stats.chisqprob(kai2(_A,E).sum(1),numpy.shape(_A)[1]-1).T
    
    rankPerson=rank(numpy.matrix(A))
    rankCluster=rank(numpy.matrix(A).T)
    if isinstance(A,list):
        rankPerson=rankPerson.tolist()[0]
        rankCluster=rankCluster.tolist()[0]

    return rankPerson,rankCluster

def dameSupportRank(A,th=0.001,maxRepeat=100,forceReturn=False):
    '''
    A:support matrix
    =>どの人が(タテ)何をどれだけ(ヨコ)Supportしているか
    
    A=[[1,0,1],
       [1,2,1],
       [0,0,1],
       [0,1,1]]
    rankPerson,rankCluster=SupportRank(A)
    
    '''    
    
    def kai2(O,E):
        '''カイ二乗検定量'''
        return numpy.multiply(O-E,O-E)/zero1(E)
    zero1=numpy.vectorize(lambda o:0.000000000000000000000001 if o==0 else o)
    floatAll=numpy.vectorize(lambda x: float(x))    
    def diffPerson(A,clusterManiacRanks):
        E=numpy.multiply(clusterManiacRanks/zero1(clusterManiacRanks.sum(1)),A.sum(1))
        return (kai2(A,E)).sum(1)         
    
    def diffCluster(A,personManiacRanks):
        E=numpy.multiply(personManiacRanks/zero1(sum(personManiacRanks)),sum(A))
        return sum(kai2(A,E))

    _A=numpy.matrix(A)
    lenPerson,lenCluster=numpy.shape(_A)

    normal=numpy.matrix([1.0/lenPerson]*lenPerson).T
    E=numpy.multiply(normal,sum(_A))
    E=numpy.multiply(E/floatAll(E.sum(1)),_A.sum(1))
    rankCluster=scipy.stats.chisqprob(sum(kai2(_A,E)),lenPerson-1)
    rankPerson=scipy.stats.chisqprob(diffPerson(_A,1-rankCluster),lenCluster-1)
    
    rankCluster/=rankCluster.sum(1)
    rankPerson/=sum(rankPerson)
    for t in xrange(maxRepeat):
        rankCluster2=scipy.stats.chisqprob(diffCluster(_A,1-rankPerson),lenPerson-1)
        rankPerson2=scipy.stats.chisqprob(diffPerson(_A,1-rankCluster2),lenCluster-1)
        rankCluster2/=rankCluster2.sum(1)
        rankPerson2/=sum(rankPerson2)
        if numpy.linalg.norm(rankCluster-rankCluster2,2)<th and numpy.linalg.norm(rankPerson-rankPerson2,2)<th:
            break
        rankCluster=rankCluster2
        rankPerson=rankPerson2
    else:
        if not forceReturn:
            raise Exception("%d回の計算では無理" % maxRepeat)
    
    rankPerson=rankPerson.T
    rankPerson=1-rankPerson
    if isinstance(A,list):
        rankPerson=rankPerson.tolist()[0]
        rankCluster=rankCluster.tolist()[0]

    return rankPerson,rankCluster
